library("here")
library("sjPlot")
library("tidyverse")
library("lme4")
library("viridis")
library("lmerTest")
library("ggplot2")
library("gridExtra")
library("gt")
library("ggthemes")

source(here("p4_analysis", "analysis_functions.R"))

p4path <- here("p4_analysis", "outputs")
p3path <- here("p3_methods", "outputs")

1. MMRR

1.1 Individual sampling

1.1.1 Summary plots

mmrr_ind <- format_mmrr(here(p3path, "mmrr_indsampling_results.csv"))

# overall error
summary_hplot(mmrr_ind, "ratio_ae", colpal = "viridis", direction = -1)

summary_hplot(mmrr_ind, "env_ae", colpal = "viridis", direction = -1)

summary_hplot(mmrr_ind, "geo_ae", colpal = "viridis", direction = -1)

# bias
summary_hplot(mmrr_ind, "ratio_err", divergent = TRUE)

summary_hplot(mmrr_ind, "comboenv_err", divergent = TRUE)

summary_hplot(mmrr_ind, "geo_err", divergent = TRUE)

1.1.2 Model summaries

run_lmer(mmrr_ind, "ratio_ae", filepath = here(p4path, "MMRR_individual_ratioerr.csv"))
[[1]]
Linear mixed effect model
statistic ~ nsamp + sampstrat + K + m + phi + H + r + (1 | seed)
Predictors Fixed Effects Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp −0.0003 8.0000 8.0000 1 15.3480K 1531.810409 1.9 × 10−319***
sampstrat 0.1078 1.5728 0.5243 3 15.3480K 100.385844 2.4 × 10−64***
K 0.0126 0.6090 0.6090 1 15.3480K 116.601181 4.4 × 10−27***
m 0.0523 10.5133 10.5133 1 15.3480K 2013.030919 0.0***
phi −0.0020 0.0155 0.0155 1 15.3480K 2.971532 0.085
H 0.0096 0.3544 0.3544 1 15.3480K 67.857541 1.9 × 10−16***
r −0.0090 0.3123 0.3123 1 15.3480K 59.800103 1.1 × 10−14***
*** p < 0.001
[[2]]
Tukey test
pairwise ~ sampstrat
Contrast Estimate SE Z ratio p
EG - G*** −0.0046 0.0016 −2.7960 0.02658997***
EG - R −0.0009 0.0016 −0.5319 0.95132229
EG - T** −0.0249 0.0016 −15.0697 0.0**
G - R 0.0037 0.0016 2.2641 0.10653142
G - T** −0.0202 0.0016 −12.2737 0.0**
R - T** −0.0240 0.0016 −14.5378 0.0**
*** p < 0.05
** p < 0.001 
run_lmer(mmrr_ind, "geo_ae", filepath = here(p4path, "MMRR_individual_geoerr.csv"))
[[1]]
Linear mixed effect model
statistic ~ nsamp + sampstrat + K + m + phi + H + r + (1 | seed)
Predictors Fixed Effects Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp −0.0001 1.2007 1.2007 1 15.3480K 1.422367e+03 9.2 × 10−298***
sampstrat 0.0525 1.3073 0.4358 3 15.3480K 5.161977e+02 1.3 × 10−319***
K 0.0060 0.1396 0.1396 1 15.3480K 1.653996e+02 1.2 × 10−37***
m 0.0215 1.7801 1.7801 1 15.3480K 2.108706e+03 0.0***
phi −0.0007 0.0018 0.0018 1 15.3480K 2.107968e+00 0.15
H 0.0056 0.1183 0.1183 1 15.3480K 1.401947e+02 3.3 × 10−32***
r 0.0000 0.0000 0.0000 1 15.3480K 5.586296e-03 0.94
*** p < 0.001
[[2]]
Tukey test
pairwise ~ sampstrat
Contrast Estimate SE Z ratio p
EG - G*** −0.0021 0.0007 −3.1344 0.009333831***
EG - R −0.0010 0.0007 −1.5145 0.428730954
EG - T** −0.0223 0.0007 −33.5784 0.0**
G - R 0.0011 0.0007 1.6199 0.367294910
G - T** −0.0202 0.0007 −30.4440 0.0**
R - T** −0.0213 0.0007 −32.0639 0.0**
*** p < 0.01
** p < 0.001 
run_lmer(mmrr_ind, "env_ae", filepath = here(p4path, "MMRR_individual_enverr.csv"))
[[1]]
Linear mixed effect model
statistic ~ nsamp + sampstrat + K + m + phi + H + r + (1 | seed)
Predictors Fixed Effects Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp −0.0001 0.7452 0.7452 1 15.3480K 1940.537011 0.0***
sampstrat 0.0350 0.0609 0.0203 3 15.3480K 52.882904 5.4 × 10−34***
K 0.0017 0.0105 0.0105 1 15.3480K 27.383651 1.7 × 10−7***
m 0.0095 0.3480 0.3480 1 15.3480K 906.190605 1.8 × 10−193***
phi −0.0003 0.0004 0.0004 1 15.3480K 1.026801 0.31
H 0.0014 0.0078 0.0078 1 15.3480K 20.317704 6.6 × 10−6***
r −0.0026 0.0262 0.0262 1 15.3480K 68.233041 1.6 × 10−16***
*** p < 0.001
[[2]]
Tukey test
pairwise ~ sampstrat
Contrast Estimate SE Z ratio p
EG - G −0.0005 0.0004 −1.2011 0.6261246
EG - R 0.0005 0.0004 1.0594 0.7142981
EG - T*** −0.0045 0.0004 −10.1643 3.7 × 10−14***
G - R 0.0010 0.0004 2.2605 0.1074219
G - T*** −0.0040 0.0004 −8.9632 2.7 × 10−14***
R - T*** −0.0050 0.0004 −11.2237 8.0 × 10−15***
*** p < 0.001

1.1.3 Megaplots

MEGAPLOT(mmrr_ind, "ratio_ae", colpal = "viridis", direction = -1)

MEGAPLOT(mmrr_ind, "ratio_err", colpal = "viridis", divergent = TRUE)

MEGAPLOT(mmrr_ind, "comboenv_err", divergent = TRUE)

MEGAPLOT(mmrr_ind, "geo_err", divergent = TRUE)

1.2 Site sampling

1.2.1 Summary plots

mmrr_site <- format_mmrr(here(p3path, "mmrr_sitesampling_results.csv"))

# overall error
summary_hplot(mmrr_site, "ratio_ae", colpal = "viridis", direction = -1)

summary_hplot(mmrr_site, "env_ae", colpal = "viridis", direction = -1)

summary_hplot(mmrr_site, "geo_ae", colpal = "viridis", direction = -1)

# bias
summary_hplot(mmrr_site, "ratio_err", divergent = TRUE)

summary_hplot(mmrr_site, "comboenv_err", divergent = TRUE)

summary_hplot(mmrr_site, "geo_err", divergent = TRUE)

1.2.2 Model summaries

run_lmer(mmrr_site, "ratio_ae", filepath = here(p4path, "MMRR_site_ratioerr.csv"))
[[1]]
Linear mixed effect model
statistic ~ nsamp + sampstrat + K + m + phi + H + r + (1 | seed)
Predictors Fixed Effects Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp −0.0054 10.9296 10.9296 1 8.6290K 774.65997483 2.4 × 10−163***
sampstrat 0.2519 1.7532 0.8766 2 8.6290K 62.13047517 1.6 × 10−27***
K 0.0003 0.0002 0.0002 1 8.6290K 0.01177688 0.910
m 0.0196 0.8264 0.8264 1 8.6290K 58.57599989 2.2 × 10−14***
phi −0.0053 0.0618 0.0618 1 8.6290K 4.37735932 0.036**
H 0.0115 0.2873 0.2873 1 8.6290K 20.36485607 6.5 × 10−6***
r −0.0184 0.7299 0.7299 1 8.6290K 51.73067018 6.9 × 10−13***
*** p < 0.001
** p < 0.05
[[2]]
Tukey test
pairwise ~ sampstrat
Contrast Estimate SE Z ratio p
EG - EQ*** −0.0340 0.0031 −10.8514 2.7 × 10−14***
EG - R** −0.0101 0.0031 −3.2163 3.7 × 10−3**
EQ - R*** 0.0239 0.0031 7.6351 8.6 × 10−14***
*** p < 0.001
** p < 0.01 
run_lmer(mmrr_site, "geo_ae", filepath = here(p4path, "MMRR_site_geoerr.csv"))
[[1]]
Linear mixed effect model
statistic ~ nsamp + sampstrat + K + m + phi + H + r + (1 | seed)
Predictors Fixed Effects Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp −0.0011 0.4225 0.4225 1 8.6290K 131.605173 3.0 × 10−30***
sampstrat 0.0702 0.1566 0.0783 2 8.6290K 24.386603 2.7 × 10−11***
K 0.0329 2.3423 2.3423 1 8.6290K 729.634765 2.5 × 10−154***
m 0.1628 57.2736 57.2736 1 8.6290K 17841.045336 0.0***
phi −0.0069 0.1030 0.1030 1 8.6290K 32.099252 1.5 × 10−8***
H 0.0093 0.1876 0.1876 1 8.6290K 58.425672 2.3 × 10−14***
r 0.0036 0.0275 0.0275 1 8.6290K 8.564055 3.4 × 10−3**
*** p < 0.001
** p < 0.01
[[2]]
Tukey test
pairwise ~ sampstrat
Contrast Estimate SE Z ratio p
EG - EQ*** 0.0098 0.0015 6.5604 1.6 × 10−10***
EG - R 0.0018 0.0015 1.2065 0.4492773
EQ - R*** −0.0080 0.0015 −5.3539 2.6 × 10−7***
*** p < 0.001 
run_lmer(mmrr_site, "env_ae", filepath = here(p4path, "MMRR_site_enverr.csv"))
[[1]]
Linear mixed effect model
statistic ~ nsamp + sampstrat + K + m + phi + H + r + (1 | seed)
Predictors Fixed Effects Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp −0.0021 1.5946 1.5946 1 8.6290K 907.15865165 1.4 × 10−189***
sampstrat 0.0968 0.2031 0.1016 2 8.6290K 57.77799916 1.2 × 10−25***
K −0.0003 0.0002 0.0002 1 8.6290K 0.09364113 0.760
m 0.0056 0.0679 0.0679 1 8.6290K 38.64283865 5.3 × 10−10***
phi −0.0018 0.0074 0.0074 1 8.6290K 4.18983195 0.041**
H 0.0027 0.0153 0.0153 1 8.6290K 8.69527846 3.2 × 10−3*
r −0.0066 0.0936 0.0936 1 8.6290K 53.23325588 3.2 × 10−13***
*** p < 0.001
** p < 0.05
* p < 0.01
[[2]]
Tukey test
pairwise ~ sampstrat
Contrast Estimate SE Z ratio p
EG - EQ*** −0.0116 0.0011 −10.4807 2.7 × 10−14***
EG - R** −0.0035 0.0011 −3.1708 4.3 × 10−3**
EQ - R*** 0.0081 0.0011 7.3099 8.2 × 10−13***
*** p < 0.001
** p < 0.01

1.2.3 Megaplots

MEGAPLOT(mmrr_site, "ratio_ae", colpal = "viridis", direction = -1)

MEGAPLOT(mmrr_site, "ratio_err", divergent = TRUE)

MEGAPLOT(mmrr_site, "comboenv_err", divergent = TRUE)

MEGAPLOT(mmrr_site, "geo_err", divergent = TRUE)

2. GDM

2.1 Individual sampling

2.1.1 Summary plots

gdm_ind <- format_gdm(here(p3path, "gdm_indsampling_results.csv"))

# overall error
summary_hplot(gdm_ind, "ratio_ae", colpal = "viridis", direction = -1)

summary_hplot(gdm_ind, "env_ae", colpal = "viridis", direction = -1)

summary_hplot(gdm_ind, "geo_ae", colpal = "viridis", direction = -1)

# bias
summary_hplot(gdm_ind, "ratio_err", divergent = TRUE)

summary_hplot(gdm_ind, "comboenv_err", divergent = TRUE)

summary_hplot(gdm_ind, "geo_err", divergent = TRUE)

2.1.2 Model summaries

run_lmer(gdm_ind, "ratio_ae", filepath = here(p4path, "GDM_individual_ratioerr.csv"))
[[1]]
Linear mixed effect model
statistic ~ nsamp + sampstrat + K + m + phi + H + r + (1 | seed)
Predictors Fixed Effects Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp −0.0003 9.4335 9.4335 1 15.3430K 8.826316e+02 1.2 × 10−188***
sampstrat 0.1099 0.7896 0.2632 3 15.3430K 2.462613e+01 6.9 × 10−16***
K 0.0231 2.0533 2.0533 1 15.3430K 1.921097e+02 2.0 × 10−43***
m 0.0799 24.4973 24.4973 1 15.3430K 2.292047e+03 0.0***
phi −0.0003 0.0002 0.0002 1 15.3430K 2.244788e-02 0.88
H 0.0020 0.0153 0.0153 1 15.3430K 1.435012e+00 0.23
r −0.0057 0.1237 0.1237 1 15.3430K 1.157134e+01 6.7 × 10−4***
*** p < 0.001
[[2]]
Tukey test
pairwise ~ sampstrat
Contrast Estimate SE Z ratio p
EG - G*** −0.0135 0.0024 −5.7327 5.9 × 10−8***
EG - R*** −0.0097 0.0024 −4.0895 2.5 × 10−4***
EG - T*** −0.0198 0.0024 −8.3745 4.1 × 10−14***
G - R 0.0039 0.0024 1.6429 0.35443857
G - T** −0.0062 0.0024 −2.6417 0.04111055**
R - T*** −0.0101 0.0024 −4.2845 1.1 × 10−4***
*** p < 0.001
** p < 0.05 
run_lmer(gdm_ind, "geo_ae", filepath = here(p4path, "GDM_individual_geoerr.csv"))
[[1]]
Linear mixed effect model
statistic ~ nsamp + sampstrat + K + m + phi + H + r + (1 | seed)
Predictors Fixed Effects Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp −0.0012 106.3637 106.3637 1 15.3430K 2702.1766883 0.0***
sampstrat 0.6274 17.5095 5.8365 3 15.3430K 148.2767520 1.0 × 10−94***
K 0.0235 2.1250 2.1250 1 15.3430K 53.9860038 2.1 × 10−13***
m −0.0852 27.8979 27.8979 1 15.3430K 708.7482021 1.1 × 10−152***
phi 0.0073 0.2059 0.2059 1 15.3430K 5.2302675 0.022**
H −0.0033 0.0413 0.0413 1 15.3430K 1.0486234 0.310
r 0.0020 0.0148 0.0148 1 15.3430K 0.3751083 0.540
*** p < 0.001
** p < 0.05
[[2]]
Tukey test
pairwise ~ sampstrat
Contrast Estimate SE Z ratio p
EG - G*** 0.0899 0.0045 19.8550 0.0***
EG - R*** 0.0204 0.0045 4.5115 3.8 × 10−5***
EG - T*** 0.0493 0.0045 10.8910 4.0 × 10−14***
G - R*** −0.0695 0.0045 −15.3432 0.0***
G - T*** −0.0406 0.0045 −8.9659 2.7 × 10−14***
R - T*** 0.0289 0.0045 6.3789 1.1 × 10−9***
*** p < 0.001 
run_lmer(gdm_ind, "env_ae", filepath = here(p4path, "GDM_individual_enverr.csv"))
[[1]]
Linear mixed effect model
statistic ~ nsamp + sampstrat + K + m + phi + H + r + (1 | seed)
Predictors Fixed Effects Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp −0.0003 5.3338 5.3338 1 15.3430K 1716.7973944 0.0***
sampstrat 0.0999 0.2357 0.0786 3 15.3430K 25.2933088 2.6 × 10−16***
K −0.0034 0.0432 0.0432 1 15.3430K 13.8946308 1.9 × 10−4***
m 0.0062 0.1486 0.1486 1 15.3430K 47.8286525 4.8 × 10−12***
phi 0.0043 0.0716 0.0716 1 15.3430K 23.0482092 1.6 × 10−6***
H 0.0009 0.0030 0.0030 1 15.3430K 0.9623282 0.33
r −0.0040 0.0619 0.0619 1 15.3430K 19.9363228 8.1 × 10−6***
*** p < 0.001
[[2]]
Tukey test
pairwise ~ sampstrat
Contrast Estimate SE Z ratio p
EG - G −0.0026 0.0013 −2.0550 0.1681013
EG - R*** −0.0046 0.0013 −3.6538 1.5 × 10−3***
EG - T** −0.0106 0.0013 −8.3556 4.3 × 10−14**
G - R −0.0020 0.0013 −1.5990 0.3791273
G - T** −0.0080 0.0013 −6.3009 1.8 × 10−9**
R - T** −0.0060 0.0013 −4.7013 1.5 × 10−5**
*** p < 0.01
** p < 0.001

2.1.3 Megaplots

MEGAPLOT(gdm_ind, "ratio_ae", colpal = "viridis", direction = -1)

MEGAPLOT(gdm_ind, "ratio_err", divergent = TRUE)

MEGAPLOT(gdm_ind, "comboenv_err", divergent = TRUE)

MEGAPLOT(gdm_ind, "geo_err", divergent = TRUE)

2.1.4 Prop NA

Confirm that the distribution of NAs is as expected and the proportions are small

MEGAPLOT(gdm_ind, "geo_coeff", aggfunc = "prop_na", colpal = "mako")

2.2 Site sampling

2.2.1 Summary plots

gdm_site <- format_gdm(here(p3path, "gdm_sitesampling_results.csv"))

# overall error
summary_hplot(gdm_site, "ratio_ae", colpal = "viridis", direction = -1)

summary_hplot(gdm_site, "env_ae", colpal = "viridis", direction = -1)

summary_hplot(gdm_site, "geo_ae", colpal = "viridis", direction = -1)

# bias
summary_hplot(gdm_site, "ratio_err", divergent = TRUE)

summary_hplot(gdm_site, "comboenv_err", divergent = TRUE)

summary_hplot(gdm_site, "geo_err", divergent = TRUE)

2.2.2 Model summaries

run_lmer(gdm_site, "ratio_ae", filepath = here(p4path, "GDM_site_ratioerr.csv"))
[[1]]
Linear mixed effect model
statistic ~ nsamp + sampstrat + K + m + phi + H + r + (1 | seed)
Predictors Fixed Effects Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp −0.0064 14.4681 14.4681 1 8.3500K 5.407432e+02 5.9 × 10−116***
sampstrat 0.3179 8.3211 4.1605 2 8.3500K 1.554993e+02 5.0 × 10−67***
K 0.0000 0.0000 0.0000 1 8.3500K 2.777680e-05 1.00
m 0.0231 1.1117 1.1117 1 8.3500K 4.155104e+01 1.2 × 10−10***
phi −0.0049 0.0503 0.0503 1 8.3500K 1.879243e+00 0.17
H 0.0007 0.0009 0.0009 1 8.3500K 3.448557e-02 0.85
r −0.0137 0.3905 0.3905 1 8.3500K 1.459672e+01 1.3 × 10−4***
*** p < 0.001
[[2]]
Tukey test
pairwise ~ sampstrat
Contrast Estimate SE Z ratio p
EG - EQ*** −0.0760 0.0044 −17.3557 0.0***
EG - R*** −0.0272 0.0044 −6.1624 2.1 × 10−9***
EQ - R*** 0.0488 0.0044 11.1986 7.3 × 10−15***
*** p < 0.001 
run_lmer(gdm_site, "geo_ae", filepath = here(p4path, "GDM_site_geoerr.csv"))
[[1]]
Linear mixed effect model
statistic ~ nsamp + sampstrat + K + m + phi + H + r + (1 | seed)
Predictors Fixed Effects Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp −0.0102 37.3428 37.3428 1 8.3501K 147.9671384 9.3 × 10−34***
sampstrat 1.1474 145.0126 72.5063 2 8.3501K 287.2989940 2.2 × 10−121***
K 0.1451 43.9976 43.9976 1 8.3500K 174.3360187 2.1 × 10−39***
m 0.6239 813.6420 813.6420 1 8.3500K 3223.9742765 0.0***
phi 0.0080 0.1332 0.1332 1 8.3500K 0.5279855 0.470
H 0.0041 0.0344 0.0344 1 8.3500K 0.1361223 0.710
r 0.0201 0.8420 0.8420 1 8.3501K 3.3362023 0.068
*** p < 0.001
[[2]]
Tukey test
pairwise ~ sampstrat
Contrast Estimate SE Z ratio p
EG - EQ*** 0.2882 0.0134 21.4290 0.0***
EG - R 0.0229 0.0136 1.6897 0.2091145
EQ - R*** −0.2652 0.0134 −19.8333 0.0***
*** p < 0.001 
run_lmer(gdm_site, "env_ae", filepath = here(p4path, "GDM_site_enverr.csv"))
[[1]]
Linear mixed effect model
statistic ~ nsamp + sampstrat + K + m + phi + H + r + (1 | seed)
Predictors Fixed Effects Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp −0.0054 10.2393 10.2393 1 8.3500K 625.37735767 3.7 × 10−133***
sampstrat 0.2722 1.5358 0.7679 2 8.3500K 46.90111271 5.6 × 10−21***
K 0.0003 0.0002 0.0002 1 8.3500K 0.01089743 0.920
m 0.0236 1.1665 1.1665 1 8.3500K 71.24432149 3.7 × 10−17***
phi 0.0010 0.0022 0.0022 1 8.3500K 0.13328216 0.720
H 0.0051 0.0544 0.0544 1 8.3500K 3.32312411 0.068
r −0.0149 0.4610 0.4610 1 8.3500K 28.15629246 1.1 × 10−7***
*** p < 0.001
[[2]]
Tukey test
pairwise ~ sampstrat
Contrast Estimate SE Z ratio p
EG - EQ*** −0.0307 0.0034 −8.9621 2.1 × 10−14***
EG - R*** −0.0267 0.0035 −7.7294 5.0 × 10−14***
EQ - R 0.0040 0.0034 1.1684 0.4721507
*** p < 0.001

2.2.3 Megaplots

MEGAPLOT(gdm_site, "ratio_ae", colpal = "viridis", direction = -1)

MEGAPLOT(gdm_site, "ratio_ae", colpal = "viridis", direction = -1)

MEGAPLOT(gdm_site, "ratio_ae", colpal = "viridis", direction = -1)

MEGAPLOT(gdm_site, "ratio_err", divergent = TRUE)

MEGAPLOT(gdm_site, "comboenv_err", divergent = TRUE)

MEGAPLOT(gdm_site, "geo_err", divergent = TRUE)

2.2.4 Prop NA

Confirm that the distribution of NAs is as expected and the proportions are small

MEGAPLOT(gdm_site, "geo_coeff", aggfunc = "prop_na", colpal = "mako")

3. Comparison of error distirbutions for MMRR and GDM

mmrr_ind %>%
  filter(sampstrat != "full") %>%
  ggplot(aes(x = ratio_ae, fill = m, colour = m)) +
  geom_density(alpha = 0.5) +
  theme_few() +
  scale_fill_viridis_d(direction = -1, end = 0.7, begin = 0.3, option = "mako") +
  scale_color_viridis_d(direction = -1, end = 0.7, begin = 0.3, option = "mako")

gdm_ind %>%
  filter(sampstrat != "full") %>%
  ggplot(aes(x = ratio_ae, fill = m, colour = m)) +
  geom_density(alpha = 0.5) +
  theme_few() +
  scale_fill_viridis_d(direction = -1, end = 0.7, begin = 0.3, option = "mako") +
  scale_color_viridis_d(direction = -1, end = 0.7, begin = 0.3, option = "mako")
## Warning: Removed 5 rows containing non-finite values (stat_density).

4. Compare results of MMRR and GDM

plot(mmrr_ind$geo_coeff, gdm_ind$geo_coeff, col = gdm_ind$m)
legend("topleft", c("0.25", "1"), col = c("black", "red"), pch = 1)

plot(mmrr_ind$comboenv_coeff, gdm_ind$comboenv_coeff)

df <- data.frame(mmrr_ind, 
                 geo_mg = gdm_ind$geo_coeff/mmrr_ind$geo_coeff, 
                 env_mg = gdm_ind$comboenv_coeff/mmrr_ind$comboenv_coeff,
                 ratio_mg = gdm_ind$ratio/mmrr_ind$ratio )
summary_hplot(df, "geo_mg")

summary_hplot(df, "env_mg")

summary_hplot(df, "ratio_mg")

MEGAPLOT(df, "geo_mg")

MEGAPLOT(df, "env_mg")

MEGAPLOT(df, "ratio_mg")

plot(mmrr_site$geo_coeff, gdm_site$geo_coeff, col = gdm_site$m)
legend("topleft", c("0.25", "1"), col = c("black", "red"), pch = 1)

plot(mmrr_site$comboenv_coeff, gdm_site$comboenv_coeff)

df <- data.frame(mmrr_site, 
                 geo_mg = gdm_site$geo_coeff/mmrr_site$geo_coeff, 
                 env_mg = gdm_site$comboenv_coeff/mmrr_site$comboenv_coeff,
                 ratio_mg = gdm_site$ratio/mmrr_site$ratio )
summary_hplot(df, "geo_mg")

summary_hplot(df, "env_mg")

summary_hplot(df, "ratio_mg")

summary_hplot(mmrr_ind, "comboenv_coeff")

summary_hplot(gdm_ind, "comboenv_coeff")

summary_hplot(mmrr_ind, "geo_coeff")

summary_hplot(gdm_ind, "geo_coeff")

Figures

arrange_figures(gdm_ind, mmrr_ind, "RAE",
                "A. GDM", "B. MMRR",
                colpal = "viridis", direction = -1, dig = 2, maxv = 0.21, minv = 0.03)

arrange_figures(gdm_site, mmrr_site, "RAE", 
                "A. GDM", "B. MMRR",
                colpal = "viridis", direction = -1, dig = 2, maxv = 0.3, minv = 0.07)